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Abstract 

We present numerical results concerning the solution of the time-harmonic Maxwell 
equations discretized by discontinuous Galerkin methods. In particular, a numeri- 
cal study of the convergence, which compares different strategies proposed in the 
literature for the elliptic Maxwell equations, is performed in the two-dimensional 
case. 

Key words: time-harmonic Maxwell's equations, discontinuous Galerkin methods. 



1 Introduction 



This work is concerned with the numerical solution of the time-harmonic 
Maxwell equations discretized by discontinuous Galerkin methods on unstruc- 
tured meshes. Our motivation for using a discontinuous Galerkin method is 
the enhanced flexibility compared to the conforming edge element method 
[12]: for instance, dealing with non- conforming meshes is straightforward and 
the choice of the local approximation space is not constrained. Nonetheless, 
before taking full advantage of these features, it is required to carefully study 
the basic ingredients of the method such as the choice of the numerical flux at 
the interface between neighboring elements. In the context of time-harmonic 
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problems, the design of efficient solution strategies for the resulting sparse 
linear systems is an equally important question. 



Previous works have shown convergence results for discontinuous Galerkin 
methods applied to the time-harmonic Maxwell equations, studied in the form 
of second-order vector wave equations. Most of these works use a mixed for- 
mulation [13,11] but discontinuous Galerkin methods on the non-mixed for- 
mulation have recently been proved to converge (interior penalty technique 
[10,1] as well as the local discontinuous Galerkin method [1]). The conver- 
gence properties of these methods in the time-domain case have been studied 
in [6] when using a centered flux and in [9] when using an upwind flux. The 
case of the upwind flux has been analyzed in [8] and [7] for the time-harmonic 
problems and the convergence has been proved only for a perturbed problem. 
The general case of Friedrichs systems and the elliptic Maxwell equations in 
particular has been treated in [4] and [5]. However, to our knowledge, no di- 
rect convergence analysis on the first-order time-harmonic system (1) has been 
conducted so far, which should be useful, for instance, when using an upwind 
flux (see subsection 2.3). The main contribution of this work is a numerical 
study of the convergence of discontinuous Galerkin methods based on centered 
and upwind fluxes applied to the first-order time-harmonic Maxwell system in 
the two-dimensional case. 

The paper is organized as follows: in section 2 we present the discretization 
method as well as the different kind of fluxes considered. In section 3, the 
convergence properties are recalled in the case of the elliptic Maxwell equations 
and the solvability of the discrete perturbed problem is analyzed in the case 
of the centered flux. In section 4 the numerical convergence is studied and 
confronted to the theoretical convergence order. Different numerical fluxes are 
then compared on two distinct examples. 



2 Discretization of the first-order time-harmonic Maxwell system 

2.1 Formulation of the continuous problem 

The system of non-dimensionalized time-harmonic Maxwell's equations can 
be written in the following form: 



where E and H are the unknown electric and magnetic fields and J is a known 
current source. The parameters e r and \i r are respectively the complex- valued 



iue r E — curl H = — J, 
\ijjjJt, r H + curl E = 0, 
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relative dielectric permittivity (integrating the electric conductivity) and the 
relative magnetic permeability; we consider here the case of linear isotropic 
media. The angular frequency of the problem is given by u. We solve Equa- 
tions (1) in a bounded domain Q, and on its boundary dfl = T a U T m , we 
impose the following boundary conditions: 

- a perfect electric conductor condition on T m , ie: n x E = on T m , 

- a Silver-Miiller (first-order absorbing boundary) condition on T a , ie: (2) 
nxE + nx(nxH)=nx E inc + n x (n x H inc ) on T a . 

The vectors E mc and H mc represent the components of an incident electro- 
magnetic wave. We can further rewrite (l) + (2), assuming J equals to 0, under 
the following form: 



iuG W + G x d x W + G y d y W + G z d z W = in Q, 
(M Fm - G n )W = on T m , 
[ (M Fa - G n )(W - W inc ) = on r.. 



(3) 



where W 



E 
H, 



is the new unknown vector and Gq 



£rh 03x3 



y0 3 x3 At 



Denot- 



ing by (e x , e y , e z ) the canonical basis of R , the matrices Gi with I e {x, y, z} 
are given by: 



3x3 N e i 

Gi = | | where for a vector n. Nn 

N** 3x3 , 



n 2 



-n, 



n z n x 



1 Tly n x 



/ 



In the following we denote by Gn the sum G x n x + G y n y + G z n z and by G^ 
and G^ its positive and negative partQ We also define |Gn| = G n — G n . 
In order to take into account the boundary conditions, the matrices Mp m and 
M-p a are given by: 



Mr 



^ 3X 3 N n 

\-N n 3x3/ 



and M Ta = \G n \. 



See [3] for further details on the derivation of this formulation. 



1 If Gn = TAT 1 is the eigenfactorization then G n = TA^T 1 where A + (resp. 
A~) only gathers the positive (resp. negative) eigenvalues. 
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2.2 Discretization 



Let Vth denote a discretization of the domain Vt into a union of conforming 
elements (tetrahedral or hexahedral elements) 

n h = u k. 

KeT h 



Hy 



of (3) in V h x V h where 



We look for the approximate solutions Wh = 
the function space 14 is defined by: 

V h = {v e[L 2 (n)] 3 /vk eT h , V\ K eV(K)}. (4) 

The term V(K) denotes a space of polynomial functions on the element K. 
We take the scalar product of the first equation of (3) by a sufficiently smooth 
vector field V and we integrate over an element K of the mesh T h : 

Jiu}(G W) t Vdx + jl ]T ddiw] Vdx = 0. 

\le{x,y,z} ) 

By using Green's formula we obtain a weak formulation involving a boundary 
term. This term is replaced in discontinuous Galerkin methods by a func- 
tion §qk which is usually referred as the numerical flux (see also Ern and 
Guermond [4,5]); the aim is then to determine W h in V h x V h such that: 

/ kjiGoWhfVdx- f Wl\ G^mdx+f (* dK (W h )) t V = 0, 

VV e Vft x V h . 

(5) 

In order to couple the element K with its neighbors for ensuring the consis- 
tency of the discretization, this numerical flux can be defined in the following 

way: 



f IfkSfIWhI + iFKGnAWh} if F e T°, 
\{M F , K + I FK Gn F )W h if F e r m , 



X -{M F , K + I FK Gn F )W h - 1 -{M F , K - I FK G nF )W™ if F e T a , 

(6) 

where T , T a and T m respectively denote the set of interior faces, the set of 
faces on r a and the set of faces on r m . I F k stands for the incidence matrix 
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between oriented faces and elements whose entries are given by: 



IFK 



if the face F does not belong to element K, 

1 if F e K and their orientations match, 

— 1 if F G K and their orientations do not match. 



We also define respectively the jump and the average of V on a face F shared 
by two elements K and K: 



IV} = I FK V K + I Fk V k and {V} = ^(V K + V k ). 

Finally, the matrix S F allows to penalize the jump of a field or of some com- 
ponents of this given field on the face F and the matrix Mp,K to be defined 
later insures the asymptotic consistency with the boundary conditions of the 
continuous problem. 



2. 3 Choice of the numerical flux 



In this study, we aim at comparing the properties of three classical numerical 
fluxes: 



- a centered flux (see [6] for the time-domain equivalent). In this case S F = 
for all the faces F and, for the boundary faces, we use: 



M, 



( 



IFK 



F,K 



03x3 N n 



NL 



3x3 



Gn F I if Fer a . 



if f e r m , 



an upwind flux (see [4,14]). In this case: 



V 



ufNnNtn 3x3 
3 x3 a^N t n N nj 



M, 



F,K 



VFN nF N t riF I FK Nn l 



\ -IfkN TIf 







3x3 



vf € r m , 



with arf, ap 7 and rjp equals to 1/2 for homogeneous media. The definition of 
M FK for F in T a is identical to the centered case. 

- a partially penalized upwind flux (local Discontinuous Galerkin method, 
see [2]). This flux is characterized by a penalization coefficient given by: 



Sp = W f^„, 0\ Mfk= VFh-/N nF N kF I FK N nF | r r 
I -IfkN^ 3x3 
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The definition of M FK for F in T a is also identical to the centered case. 



3 Convergence properties of the discretized problem 

We are interested in assessing these numerical fluxes for the discretization of 
(3). Firstly, we want the best asymptotic convergence order in L 2 -norm for 
the electric and magnetic field for a fixed polynomial order approximation 
on an unstructured mesh. Secondly, a minimal numerical dispersion is also 
needed. In the following we will focus on the first criterion. The asymptotic 
convergence order in L 2 -norm between the exact solution (E, H) and the 
approximate solution (E h) H h ) corresponds to the largest real coefficients (3 
and 7 such that: 

3C u C 2 ,ho > 0, V/i > ho, \\E-E h \\ L 2 (n) < dtf and \\H-H h \\ L2(n) < C 2 h\ 

(7) 

where h is the mesh size. Let us note that in the numerical examples proposed 
in Section 4, we will often equivalently consider the evolution of the norm of 
the error against the square root of the number of degrees of freedom (dofs), 
in order to deduce coefficients (3 and 7. 

We first recall in Table 1 below the theoretical convergence order for the elliptic 
Maxwell equations [4,5], for a sufficiently smooth solution and when the local 
function space V{K) is [Pk(K)] 3 i.e. the space of vectors whose components 
are polynomials of order at most k. When using the flux with a penalization 
of E, similar convergence results are proved for the time-harmonic Maxwell 
equations in [1]. 



flux 


centered 


upwind 


penalization of E 


field E 


k 


fc + 1/2 


k + 1 


field H 


k 


fc + 1/2 


k 



Table f 

Theoretical convergence order for the elliptic Maxwell equations. 



3.1 Solution of the discretized perturbed problem 

A few comments need to be stated concerning the convergence properties of 
such a scheme applied to the first-order formulation of the time-harmonic 
Maxwell equations. First of all, the case of the upwind flux has been analyzed 
in [8] for the perturbed Maxwell problem, that is when iuj is replaced by 
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v + iu with v a strictly positive parameter. For a sufficiently regular solution 
the norm of the error behaves as h p+1 ^ 2 where h is the mesh parameter. 



The case of the centered flux has been studied in [6] for the time-domain 
Maxwell equations and in this case the norm of the error behaves as hP where 
h is the mesh parameter. For the time-harmonic equations no convergence 
proofs are available so far. We can only study here the solvability of the discrete 
problem in the case of a perturbed problem (we replace iuj by iu+u with v > 0) 
following an idea used by Helluy [7] in the case of the upwind flux. In the case 
of the perturbed problem and assuming homogeneous boundary conditions, 
the formulation can be simply written as: 



J Find Wh in Vh x Vh such that: 

\ a(W h , V) + b(W h , V) = 0, VV e V h x V h , 



with, VC7, V G V h x V h : 



a(U,V)=[ ((iu; + u)G UyVdv+ £ / ( l -\G nF \U) Vds 
+ E / (\M F , K U)Vds+ J2 I {SrWlfWiFds, 



Fer™ r x 7 Fer ' 



and: 



t 

Vdv 



b(U,V)= E / f E GA(U) 

K£T h JK \le{x,y,z} j 

- E / (\lFKGn F U)Vds ( 10 ) 

FGr a ur m 

- E [ (GnAUjf {V}ds. 
We have the following result: 

Proposition 1 ITie solution of problem (8) zs nu//. 

Proof First, considering the fact that the matrices |Gn F |, Sp, 3?(C ) and 
— ^(Gq) are hermitian and denoting by Ti.{M FK ) the hermitian part of M F>K 
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for F in T m , which is equal to 



VfNu f NL 3x3 



J3x3 



! one has: 

3 x3,' 



&(a(W h , W h )) = / ((^(G ) - u%(G ))W h y W h dv 



+ 



Fer° JF 

+ E / {\\Gn F \wSwids 

Fera JF\2 J 

+ E l F {^n{M F , K )w^)WHds. 



Then, we rewrite using the corresponding Green identity an equivalent expres- 
sion of the sesquilinear form b: 



b(U,V) = - £ 

K€T h 



j K Ut [ E GA(V)\dv- 

\le{x,y,z} ) 

E / (lFKGn F U lK yV lK ds 



FedK 



^ jf (GnAUjf {V}ds, VC7, V e 14 x 14. 



Fer° 



By noticing that on a face F e T° separating two elements K and K: 

(GnAUmVj + (GnjC/])*{V} = (I FK Gn F U\ K yV\ K 

+ (I F K G n F U^yV^, 

which is in part due to the fact that Gn F is hermitian, one deduces: 



b(U, V) 



KeT h 

+ E 



E / ^ ( E ^(v) ] dv 

cer h Jk \le{x,y,z} J 



Fer a 



.\l F KGn F U\V(h 

F \Z 



+ E / F (Gn f {[/}) ( [V]ds, VC7, V e V h x 14. 



(12) 



(13) 



Fer° 



Thus, it is now straightforward to see that 6 is anti-hermitian and conse- 



quently: 




®(a(W h , W h ) + b(W h , W h )) = f ((i/ft(G ) - w9(G )) W h dv 



From (8), &(a(W h , W h ) + 6(W h , W h )) is also equal to zero. As i/ft(G ) - 
u;$5(G ) is positive definite and |Gn F |, Si? an d TL(M F>K ) are positive, the 



4 Numerical results 

In the first part of this section we will present a numerical comparison of 
different fluxes for a very simple test case and different kind of meshes. In 
the second part, the results on a less trivial problem are compared to those 
obtained with the plane wave example. 

We consider the case of an electric transverse wave in the plane (0,x,y). In 
this case the components E z , H x and H y are zero. We numerically simulate 
the propagation of a plane wave in vacuum where the incident wave is given 
by (E™ c , Ey lc , H" 1C ) = exp(— iux)(Q, 1, 1). The computational domain is the 
unit square Q =]0; 1[ 2 and a Silver- Miiller boundary condition is imposed on 
the whole boundary, that is T a = dfl and T rn = 0. The parameters e r and 
fi r are set to 1 everywhere and we choose uj = 2ti. We numerically estimate 
the asymptotic convergence order of discontinuous Galerkin methods for the 
above problem using two different sequences of triangular meshes: 

- uniformly refined meshes. The first mesh of Figure 1(a) is uniformly 
refined resulting in the meshes of Figures 1(b) and 1(c). 

- independent meshes. We use four unstructured (quasi-uniform) indepen- 
dent meshes with an imposed maximal mesh size h (see Figure 2 for the first 
three meshes). These meshes are denoted by % for % — 1, . . . , 4 with h in a 
decreasing order. Thus T i+i is not a refinement of %. 

Our implementation of high order discontinuous Galerkin methods makes use 
of nodal basis functions with equi-spaced nodes. 



vector field Wh is null. 
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(a)/t"=f/8." ' ' (b) h = 1/16°" ' ° (c) /i = 1/32." 

Fig. 2. First three independent unstructured meshes. 

4-1 Convergence behavior using meshes obtained by uniform refinement 



Centered flux. Numerical convergence results in a logarithmic scale are 
shown on Figure 3. They clearly demonstrate the interest of higher order poly- 
nomial approximations which allow a considerable reduction of the number of 
degrees of freedom to reach the same accuracy. Table 2 summarizes numerical 
estimates (using a linear regression method) of the asymptotic convergence 
order. 




10 1 10 2 10' 10 ! 



VNumber of dofs ^Number of dofs 

(a) ||JjT— i3"^||i2 against the square root (b) \\E — Eh\\ L 2 against the square root 
of the number of dofs. of the number of dofs. 

Fig. 3. Convergence results using a centered flux. Solid lines show the evolution for 
the whole of the numerical results and dotted lines show the asymptotic tendency, 
using coefficients [3 or 7 from inequalities (7) estimated by a linear regression. 

The method based on a Pq approximation (i.e. the standard cell centered 
finite volume method) is special: the convergence order is optimal for both 
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Pa 


Pi 


P2 


^3 


E 


1.0 


1.0 


2.0 


3.0 


H 


1.0 


2.0 


3.0 


3.6 



Table 2 

Numerical convergence order using a centered flux. 

fields E and H, that is, equal to k + 1. This could be the consequence of using 
uniformly refined meshes, since a somewhat different behavior is obtained for 
independent meshes with decreasing mesh size (see subsection 4.2). For the 
other polynomial degrees, we get exactly the predicted theoretical convergence 
order in the elliptic case for E, whereas for H, this convergence order is 
optimal. Therefore, in this example, the magnetic field is better approximated 
than the electric field, when using the centered flux. 

Upwind flux. We used here the parameters ctp = ctp — r\p — 1 for each face 
F. Numerical convergence results are shown on Figure 4. Similar conclusions 
can be derived as in the centered case except that the convergence properties of 
the methods based on P and Pi interpolations are this time clearly different 
with respect to the centered case. The asymptotic convergence orders (see 
Table 3) are similar for both fields and correspond to the theory for the elliptic 
Maxwell equations. The convergence is optimal except for the case Pq, but 
nevertheless we are still above the theoretical estimates. 




10 1 10 2 10 1 10 2 



VNumbcr of dofs VNumbcr of dofs 

(a) \\H — Hh\\ L 2 against the square root (b) \\E — Eh\\ L 2 against the square root 
of the number of dofs. of the number of dofs. 

Fig. 4. Convergence results using an upwind flux. Solid lines show the evolution for 
the whole of the numerical results and dotted lines show the asymptotic tendency, 
using coefficients j3 or 7 from inequalities (7) estimated by a linear regression. 





Pa 


Pi 


P2 


P3 


E 


0.9 


1.9 


3.0 


3.9 


H 


0.9 


1.9 


3.0 


3.9 



Table 3 

Numerical convergence order using an upwind flux. 
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Penalized flux on E. We set t> = rjp = 1 for each face F. Results are shown 
on Figure 5. Table 4 summarizes the numerical estimates of the asymptotic 
convergence order. Besides the expected lack of convergence in the case Pq, 
we can notice for all the other cases ((Pk)k>o) a complementary behavior with 
respect to the centered flux, since this time we get an optimal convergence 
rate for E, but not for H . 




Hf_ 10 2 

VNumbcr of dofs \/Number of dofs 



(a) 1 1 JET — ffftH^a against the square root (b) \\E — Eh\\ L 2 against the square root 
of the number of dofs. of the number of dofs. 

Fig. 5. Convergence results using a penalized flux on E. Solid lines show the 
evolution for the whole of the numerical results and dotted lines show the asymp- 
totic tendency, using coefficients (3 or 7 from inequalities (7) estimated by a linear 
regression. 





Pa 


Pi 


P2 


^3 


E 


X 


2.0 


3.1 


3.9 


H 


X 


1.0 


2.0 


2.9 



Table 4 

Numerical convergence order using a penalized flux on E. 

4-2 Convergence behavior using independent meshes 

On Figure 6, we compare the evolution of the L 2 -norm of the error with 
the mesh size h by using the meshes (^)j=i i ...,4, for both a centered flux and 
an upwind flux, Figure 6(b) corresponds to the error for the field E while 
Figure 6(a) corresponds to the error for the field H. The results for the upwind 
flux are the same as for the uniformly refined meshes. For the centered flux, 
note the lack of convergence for the case Pq. For all the other cases the results 
remain the same as for the uniformly refined meshes. 

It is already known for time-domain problems that the centered flux combined 
to a leap-frog time integration scheme results in a non-dissipative discontin- 
uous Galerkin method (a mandatory feature for long time computations, see 
[6]). As far as time-harmonic problems are concerned, the previous results show 
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10" 



n 10 




H — I* Po centered 
h — f Pq upwind 

h — ^ Pi centered 
* — * P] upwind 

* — ► P2 centered 
P% upwind 



10° 



H io 3 




H t Po centered 

H — ^ Po upwind 

H 1 h «M 

^ ^ Pi centered 

< — > Pi upwind 

^ i Pt centered 

P2 upwind 
+ ■ ■* ft 1 



10" 



10" 



Mesh size h 



Mesh size h 



(a) ll-H" — Hh\\ L 2 against the mesh size h. (b) \\E — Eh\\ L 2 against the mesh size h. 

Fig. 6. Comparison of the convergence results between centered flux and upwind 
flux. 

that the upwind flux has better convergence properties. Nevertheless, the cen- 
tered flux remains less expensive both for time-domain and time-harmonic 
problems (arithmetic operations and memory requirements). 



4-3 Numerical comparisons on a less trivial problem 



The domain is the square [— 1; l] 2 where we have suppressed a part by inserting 
a point of coordinates (0.1,0) at it is shown on Figure 7. The properties e r 
and ii r are still homogeneous and equal to one. Appropriate non-homogeneous 
Dirichlet boundary conditions are enforced on the boundary of the domain in 
order to obtain E = (sin(27n/), sin(27nr))' as the solution. 

The mesh is not fully homogeneous as it is shown on Figure 7; it is slightly 
denser next to the point of coordinates (0.1, 0). Independent meshes have been 
used as in Subsection (4.2). 




(a) First mesh. /i max = (b) Second mesh. /i max = (c) Third mesh. /i max = 
0.32 0.16 0.32 

Fig. 7. Three first meshes used for the second example. 

The results shown on Figure 8 are in a full agreement with those obtained in 
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the case of the plane wave and independent meshes at Subsection 4.2. 




VNumber of dofs 




vNumber of dofs 



(b) Evolution of the L 2 -norm of the 
error for the H field. 



(a) Evolution of the L 2 -norm of the 
error for the E field. 

Fig. 8. Evolution of the L 2 -norm of the error against the square root of the number 
of degrees of freedom (dofs). 
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